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Computer  programs  for  high-current  beam  transport  in 
accelerators 

'  r'r 

Brer.dari  B.  Godfrey 

Mission  Research  Corporation,  Albuquerque,  NM  87106,  LSA 


Numerical  techniques  exist  for  modeling  particle  beam  dynamics  in  high  current 

acceierators  at  several  levels  of  fidelity.  Equilibria  can  be  determined  with  beam  _ 

envelope  codes,  particle  ray-tracing,  V'lasov  equilibrium  solvers,  and  single-disk  par-  r 
tide  cooes.  The  linear  stability  of  these  equilibria  is  addressed  with  dispersion  re-  y 
lations  solved  numerically,  linearized  PIC  codes,  and  single-  and  multi-component  t 

beam  centroid  programs.  Beam  nonlinear  dynamics  are  investigated  with  multidi¬ 
mensional  PIC  codes  employing  either  the  complete  electromagnetic  field  equations  _ 
or  various  approximations  to  them.  Each  of  these  options  is  discussed,  and  several_ 
examples  are  provided.  .  t:  i,  ;  i  ;■/ 

1.  INTRODUCTION  :  t 

A.  Accelerator  Types 

Historically,  high  energy  particle  accelerators  were  developed  as  research  tools  for  nuclear  and 
high  energy  physics  research.  These  accelerators  were  low  current  devices  for  which  the  beam 
dynamics  could  be  modeled  computationally  with  a  small  number  of  noninteracting  particles. 
Spacecharge,  wakefield,  and  other  corrections  were  added  to  the  calculations  as  the  currents  of 
interest  gradually  increased.  Several  sophisticated  yet  reasonably  fast-running  computer  programs 
now  exist  for  modest  current,  high  energy  accelerators;  they  are  described  in  the  chapter  by 
R.  Cooper. 

Since  the  I960's  high  current,  moderate  energy,  pulsed-diode  electron  accelerators  have  been 
buiit.  principally  for  materials  research  and  flash  radiography.  More  recently,  similar  techniques 
have  been  used  to  create  intense  light  ion  beams  for  fusion  purposes.  Self-consistently  treating 
beam-generated  fields  is  critical  for  accurate  numerical  simulation  of  high  current  diodes,  and  even 
the  earliest  steady-state  codes  included  first-principles  calculations  of  electrostatic  and  azimuthal 
magnetic  fields.  Time-dependent  diode  codes,  borrowed  from  plasma  physics,  also  include  self¬ 
fields  in  an  integral  manner.  The  chapter  by  J.  Quintenz  and  D.  Seidel  provides  details. 

High  current,  high  energy  induction  accelerators  first  were  studied  in  the  late  1960’s  as  injectors 
"or  electron  ring  accelerators.^  (Research  on  betatrons,  of  course,  goes  back  much  further.'"' 
•More  recently,  renewed  investigations  have  been  motivated  by  heavy  ion  fusion,  free  electron 
•  asers.  and  other  applications.  Treating  induction  accelerators  numerically  is  challenging  due  to 
'he  large  range  of  temporal  and  spatial  scales  involved.  Simultaneously,  self-field  effects  must  be 
.nciuded.  especially  for  low  energy  sections.  A  wide  variety  of  computational  tools,  from  simple 
-mveiope  models  to  multidimensional  particle-in-cell  (PIC)  simulations,  have  been  developed  and 
are  reviewed  in  this  chapter. 


Several  collective-effect  concepts  have  been  proposed  to  utilize  the  very  large  fields  which  can. 
in  principle,  exist  in  plasmas  to  accelerate  electrons  or  ions.  Collective  ion  acceleration.*  ^  laser 
heatwave  acceleration,®’^  and  plasma  Wakefield  acceleration^  are  examples.  Typically,  PIC  codes 
are  used  to  explore  these  ideas  numerically. 

B.  High  Power  Induction  Accelerators 

linear  induction  accelerator  achieves  high  particle  energies  by  applying  moderate  voltages 
at  successive  acceleration  gaps  only  while  the  beam  is  passing.  The  gaps  are  driven  by  ferrite 
cores  or  dielectric  lines,  which  are  in  turn  energized  by  conventional  pulse  power  sources.  Figure  1 
illustrates  schematically  a  portion  of  the  beam  line  of  a  ferrite- loaded  accelerator.  Pavlovskii  and 


FIG.  1.  Simplified  representation  of  induction  accelerator  beam  line,  showing  two  ferrite  core 
modules.  Solenoidal  magnetic  field  coils  typically  are  placed  between  gaps. 

co-workers  published  extensively  on  linear  induction  accelerators  during  the  1970’s.®  Two  half- 
terawatt  electron  accelerators  are  in  operation,  the  10  kA,  42  MeV,  ferrite  core  ATA  at  Lawrence 
Livermore  National  Laboratory^®  and  the  30  kA.  16  MeV  dielectric  line  RADLAC  at  Sandia 
-National  Laboratories. “ 

The  acceleration  gaps  can  perturb  the  beam  in  several  ways.  Radial  electric  fringe  fields 
may  cause  radial  oscillations,^^  while  the  interrupted  magnetic  image  current  may  trigger  trans¬ 
verse  displacements.^®  The  beam  breakup  instability,  associated  with  T\finQ  gap  modes,  produces 
particularly  dangerous  transverse  oscillations.^*  *®  Even  in  the  absence  of  gaps,  the  diocotron  in¬ 
stability  is  disruptive  for  low  energy  annular  beams.*®  The  resistive  wall  instability  is  unimportant, 
however.*^ 


\'arious  recirculation  schemes  are  being  considered  to  reduce  the  size,  weight,  and  cost  o: 
induction  accelerators. Figure  2  depicts  a  basic  design.  An  electron  beam  is  injected  through 


FIG.  2.  Schematic  representation  of  a  racetrack  induction  accelerator. 

an  entrance  port,  passes  several  times  through  the  acceleration  modules,  and  is  extracted  through 
the  exit  port.  Because  injection  and  extraction  with  minimal  beam  loss  is  difficult  even  concep¬ 
tually,  some  designs  employ  an  open  helix  rather  them  a  closed  racetrack  geometry.^®  Bending 
magnetic  fields  in  the  curved  sections  must  be  ramped  up  rapidly  to  keep  the  beam  in  the  drift 
tube.  Periodic  strong-focusing  fields  in  the  bends  is  sometimes  proposed  to  incretise  the  energy 
bandwidth  there. The  betatron  also  is  a  recirculating  induction  accelerator,  although  without 
discrete  gaps.  A  slowly  increasing  vertical  magnetic  both  guides  and  accelerates  the  beam.*^  In 
addition,  a  toroidal  magnetic  field,  perhaps  augmented  by  periodic  strong  focusing,  is  needed  to 
confine  high  currents  at  low  to  moderate  energies.”  Injection  and  extraction  remain  problems. 

Beam  line  curvature  introduces  the  negative  mass  instability,^®'*^  which  can  destroy  a  recircu¬ 
lating  beam  on  the  ^sec  time  scale.  Figure  3  shows  the  typical  saturated  state  of  a  long  wavelength 
negative  mass  mode  in  a  high  current  betatron.*®  The  instability  arises  because  higher  energy 
beam  particles  in  a  bend  move  to  larger  radii  and  consequently  fall  back  in  the  beam  relative  to 
lower  energy  particles.  Particle  clumps  grow  spontaneously,  as  if  like  particles  attracted.  Periodic 
strong  focusing,  mentioned  above,  slows  the  negative  mass  instability  but  triggers  a  parametric 
three-wave  instability  similar  to  the  Raman  instability  of  the  free  electron  laser.*®  Its  growth 
rate  can  be  very  large.  Single  particle  resonances  also  are  potentially  serious  in  recirculating 
accelerators,  so  special  care  must  be  taken  to  avoid  magnetic  field  errors. 


FIG.  3.  The  (a)  z-r  and  (b)  r-d  cross  sections  of  a  10  kA,  5  MeV  electron  beam  in  a  1  m  radius 
betatron  260  nsec  after  the  onset  of  a  long  wavelength  negative  mass  instability;  Bj  = 
1  kG. 


Ion-focused  transport  (IFT)  in  induction  electron  accelerators,  either  separately  or  in  con¬ 
junction  with  magnetic  transport,  has  certain  desirable  features^^  magnetic  focusing  elements 
are  unnecessary,  the  beam  breakup  instability  is  suppressed,®^  and  the  energy  bandwidth  of  mag¬ 
netic  bending  elements  is  increased.  Injection  and  extraction  for  racetrack  geometries  remain 
problems,  although  more  possible  solutions  exist.®®  Low  density  ionized  channels  can  be  created 
by  lasers,  low  energy  electron  beams,  or  electrical  discharges.  In  an  experiment  at  Sandia  National 
Laboratories,  a  10  m  ion  channel  was  produced  with  a  1  A,  500  V  electron  beam,  and  a  15  kA, 
1  MeV  electron  beam  guided  by  the  channel  through  three  90°  bends. ®^  The  most  serious  poten¬ 
tial  difficulty  for  IFT  is  the  ion  resonance  instability.®®  ®®  It  can,  at  a  minimum,  be  postponed  by 
using  heavy  ions,  like  Xenon. 

Experiments  on  the  PULSELAC  device  have  demonstrated  that  electron  neutralization  in  ion 
accelerators  also  is  effective.®^  Extensive  reviews  of  high  current  induction  accelerators  can  be 
found  in  the  proceedings  of  a  recent  NATO  Advanced  Study  Institute.®* 


C.  Numerical  Models 


Beam  dynamics  in  induction  and  other  high-current  accelerators  can  be  treated  numerically 
in  a  variety  of  approximations.  Beam  equilibria  often  are  computed  economically  with  envelope 
codes.  Steady-state  particle  ray-tracing  codes  are  useful  in  modeling  diodes  and  short  trans¬ 
port  sections.  If  the  beam  particle  equilibrium  distribution  is  approximately  known  (in  terms 
of  constants  of  motion),  then  the  beam  current  and  field  radial  profiles  are  easily  determined 
numerically.  Finally,  radially-resolved  single  disk  'C  codes  axe  used  to  obtain  the  evolution  of 
the  beam  equilibrium  when  beam  properties  vary  slowly  through  the  pulse. 

Beam  stability  properties  in  an  accelerator  can  be  determined  simply  by  solving  numerically 
an  analytical  linear  dispersion  relation,  if  it  is  known,  .\lternatively,  linearized  PIC  codes  are 


fairly  efficient  at  providing  the  growth  rates  and  field  structure  of  ’hree-dimensional  instabilities 
for  specified  one-dimensional  (typically,  radially-resolved)  equilibria.  These  two  procedures  yield 
instability  growth  rates  at  fixed  points  in  the  accelerator.  The  linear  development  of  transverse 
oscillations  due  to,  for  instance,  the  beam  breakup  instability  as  the  beam  travels  through  the 
entire  accelerator  can  be  modeled  efficiently  by  treating  the  beam  as  a  string  of  rigid  disks,  and 
the  gap  fields  as  lumped  circuit  elements.  (Early  stages  of  the  negative  mass  instability  have 
been  handled  in  the  same  way  as  well.)  Beam  thermal  effects  sometimes  can  be  introduced  by 
representing  the  beam  by  multiple  sets  of  interlaced  disks.  If  required,  the  disk  radii  can  evolve 
according  to  the  envelope  equations  mentioned  above. 

.\t  some  point,  adding  more  and  more  physics  to  simple  numerical  models  becomes  unproduc¬ 
tive,  and  PIC  codes  should  be  used  instead.  This  is  particularly  true,  if  the  nonlinear  evolution 
of  instabilities  is  important.  Though  expensive  to  employ,  PIC  codes  describe  the  interaction 
of  beam  particles  with  applied  and  self-consistent  electromagnetic  fields  in  microscopic  detail. 
Typically,  PIC  accelerator  simulations  are  performed  in  two-dimensional  r-z  geometry,  although 
three-dimensional  calculations  sometimes  are  necessary.  Single  disk  PIC  codes  resolved  in  v-Q 
can  be  used  to  investigate  long  wavelength  transverse  instabilities.  When  electromagnetic  waves 
are  unimportant,  employing  the  Darwin  field  approximation  saves  some  time.  The  frozen  field 
approximation  could  be  applied  to  plasma  wakefield  calculations. 


II.  EQUILIBRIUM  AND  QUASI-EQUILIBRIUM  COMPUTATIONS 

A.  Beam  Envelope  Codes 

In  most  cases  1  for  relativistic  electron  beams  in  accelerators,  although  not  necessarily 

in  diodes.  (V'  is  the  beam  current  normalized  to  17  kA,  7  =  (1  -  is  the  particle  relativistic 

energy  normalized  to  its  rest  mass,  and  V  is  the  particle  velocity  normalized  to  the  speed  of  light.) 

.\s  a  result  V'_  <g;  and  the  paraxial  approximation  that  changes  in  do  not  cause  changes  in 
r,  can  be  invoked,  .\pplying  the  virial  theorem  to  the  beam  transverse  dynairnics  then  yields  a 
simple  differential  equation  for  the  beam  rms  radius,  R.^® 


d  d  ^  r  -  ~iU  ^ 

dt  Ut  -rR^  47  R 


.\uxiiiary  quantities  are  the  normalized  emittance  e. 


and  the  self-field  function  U, 


U  =  (E,  -  VB»).  V'yr^  «  V- 


1  - 


1-/ 


(1) 


(3) 


Pd  is  the  particle  angular  momentum  normalized  to  me,  =  qB^lmc  is  its  gyrofrequency, 
is  any  applied  radial  acceleration,  and  /  is  the  charge  neutralization  fraction  provided  by  a 


background  piasma.  Current  neutralization  could  be  incorporated  in  Eq.  (3)  as  well.  Note  that 
Eq.  ;^1)  reduces  to  the  single  particle  radial  equation  of  motion  for  e  small. 

Equation  (1)  as  it  stands  is  dissipationless.  and  the  normalized  emittance  is  assumed  constant. 
In  reality,  single  particle  scattering  causes  a  slow  growth  in  emittance. as  do  weak  instabilities. 
In  addition,  phase-mix  damping  due  to  a  spread  in  single  particle  oscillation  frequencies  leads  to 
a  gradual  decrease  in  d‘ R  dt'  with  a  corresponding  increase  in  £.■“  These  corrections  usually  can 
be  ignored  in  accelerator  calculations. 

One  of  many  problems  to  which  Eq.  (1)  hzis  been  applied  is  the  chopping  of  the  initially  contin¬ 
uous  beam  in  the  PHERMEX  standing-wave  RF  accelerator  at  Los  Alamos  National  Laboratory.^* 
The  effects  of  apertures  and  RF  focusing  fields  as  well  as  RF  accelerating  fields  are  included.  De¬ 
spite  the  simplicity  of  the  model,  agreement  between  code  predictions  and  experimental  measure¬ 
ments  is  good.  Figure  4  illustrates  typical  output  for  an  individual  micropulse  out  of  a  sequence 
of  several.  Calculations  take  less  than  a  minute  on  a  minicomputer. 
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FIG,  4.  PHERMEX  beam  current  (/(,),  energy  (7),  radius  (rj),  and  divergence  (qd)  versus  RF 
phase  as  predicted  by  XLR8R  code. 

The  envelope  model  also  can  be  used  to  represent  one  particle  species  in  a  multispecies  hybrid 
model.  For  instance,  the  SEE  code  at  Lawrence  Livermore  National  Laboratory,  which  is  used  to 
investigate  secondary  electron  expulsion  from  an  IFT  channel,  describes  the  electron  beam  with 
an  envelope  equation  and  the  channel  electrons  with  discrete  particles. 


B.  Steady-State  Ray-Tracing  Codes 

Long-puise  electron  and  ion  diodes  are  treated  efficiently  by  ray-tracing  codes  with  applied 
and  seif-consistent  fields.^  These  codes  aiso  can  be  used  for  steadv-state  beam  transoort  in 


P.'F?  f»J  i  Wl*^  ".■  fV»  ■j'  '"Ji  ■  J  "Ji  -J  "J  "J  •Ui  ■<  ■•»  ■”»  ■  «  -j  ’■■r  ■  W  -  W  -J.- 

accelerator  segments.  At  each  iteration  a  few  hundred  beam  particles  are  transported  aiong  their 
entire  trajectories  on  the  computational  mesh  using  fields  determined  on  the  previous  iteration, 
and  their  charge  and  current  contributions  accumulated.  The  static  beam  fields  then  are  recom¬ 
puted  from  the  charge  and  current  distributions.  This  process  is  repeated  a  few  dozen  times  until 
convergence  is  achieved.  Only  about  a  minute  of  CRAY-1  computer  time  is  needed.  Often  the 
output  is  used  to  initialize  other  computer  programs.'*'’  Particle  ray-tracing  codes  are  described 
more  fully  in  the  chapter  by  J.  Quintenz  and  D.  Seidel. 

C.  Kinetic  Equilibrium  Codes 

.\ny  cylindrically  symmetric  beam  equilibrium  can  be  expressed  as  a  distribution  function 
FiH,  P,,  Pg)  of  the  particle  energy  and  canonical  axial  and  angular  momenta.*' 


//■  =  -»  ^  C>(r)  (4) 

P,  =p, -.4,(r)  ,5) 

Pi  =  r[pf  -r  As(r)]  (6) 

.\ll  the  beam  properties  can  be  determined  from  F.  The  scaler  and  vector  potentials  appear¬ 
ing  in  Eqs.  (-ij-je)  are  given  in  terms  of  the  beam  charge  and  current  densities  by  the  usual 
expressions. 


1  d  d 

-  r~o  =  -47rp  7) 

r  dr  ar 


I  d  d 

-—r—-A,  =  -47rJj  (S'l 

r  dr  dr 


d  1  d 

r - -rAi  =  -AnJ)  {9l 

dr  r  dr 

The  charge  and  current  densities  are.  in  turn,  given  by  integrals  over  the  distribution  function. 


P  =  J 

1  Fd^p 

(10) 

J,  =  l 

'  VrFd^p 

(11) 

ViFd^p 

(12) 

Lsually.  it  is  more  convenient  to  perform  the  integrals  over  the  constants  of  motion  instead  of  the 
kinetic  momenta.  The  Jacobian  of  the  transformation  is  7  rp,.  Care  must  be  taken  in  evaluating 
the  integrals  near  radial  turning  points,  where  p^  =  0. 


^  •  J»  J*  •' 

-*x  .*»  /r  i 


In  the  absence  of  instabilities,  reasonable  guesses  can  be  made  for  the  particle  distribution 
function  in  an  induction  accelerator.  Since  canonical  angular  momentum  is  conserved,  the  Pi 
distribution  is  determined  by  injector  properties.  For  a  shielded  injector  (i.e.,  Bl  =  0).  Pi  -  Q 
and  F  is  proportional  to  6(Pi].  More  generally,  the  P^  distribution  is  estimated  by  solving 


'.Aj(rj  =  P^  for  r  at  the  cathode  and  then  substituting  for  r  in  n'{r)’rB\  at  the  cathode.  The 
H  distribution  is  obtained  with  equal  ease:  H  is  approximately  constant  for  all  particles  at  the 
.njector  and  increases  by  a  constant  amount  at  each  gap.  Finally,  the  spread  in  P,  is  determined  by 
the  angular  scatter  in  the  beam  particles,  which  is  small.  One  plausible  expression  is  expjP^  T), 
where  T  =  po  ■  p\  Pi  -  is  proportional  to  the  Bennett  temperature.^*'*®  {T  can  be  made  a 
function  of  Pi.  if  a  strong  correlation  exists  between  the  particle  scatter  angle  and  radius.)  In  all. 
we  have  an  expression  of  the  sort 


F{H.  P..  Pi)  ^C6(H  -  Ho)  exp  (P.;T)  {n/rB,}  (>.4,;"'  (P,))  ( 13) 

Equations  |4i-(12)  are  easily  solved  numerically.  Density  profiles  obtained  from  the  computer 
program  ORBIT  for  four  simple  beam  distribution  functions  are  shown  in  Fig.  5.^®  Convergence 
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FIG,  3.  Sample  density  profiles  for  monoenergetic  beams,  obtained  by  solving  Eqs.  (4)-(12i  with 
the  indicated  distribution  functions. 

s  achieved  m  about  ten  iterations,  which  can  be  done  in  a  few  seconds  on  a  CRAY-1  computer. 
P.quiiibria  obtained  in  this  way  have  been  used  to  initialize  some  of  the  linear  and  nonlinear  PIC 


codes  discussed  later.  Incidentally,  the  distribution  C  ■  6{H  -  //o)  -explP^  T]  is  a  nonparaxiai 
generalization  of  the  Bennett  distribution. 

ORBIT  has  been  employed  in  magnetically  Insulated  ion  diodes  as  well.  Finite  temperature 
effects,  gcis  prefill,  and  electron  injection  with  nonzero  energy  were  investigated,  and  both  axial 
and  azimuthal  applied  magnetic  fields  were  considered. Using  ORBIT  eliminated  the  need  for 
many  of  the  assumptions  and  approximations  necessary  in  the  analytic  treatment  of  magnetic 
insulation. 

D.  Single  Disk  Particle  Codes 

The  particle  constants  of  motion  are  changed,  if  an  abrupt  transition  occurs  in  the  external 
focusing  fields;  a  dynamic  calculation  is  necessary  to  obtain  the  new  beam  equilibrium.  This 
is  particularly  simple  in  the  case  of  relativistic  beams,  for  which  the  effects  of  the  transverse 
self-fields  (e.g.,  Er  and  Bg)  nearly  cancel.  The  beam  then  is  represented  as  a  single  disk  of  a  few 
hundred  noninteracting  particles  subject  only  to  the  applied  fields.  Such  codes  can  be  very  fast, 
even  for  transport  over  long  distances.  (Single  disk  particle  codes  in  which  transverse  seif-fields 
are  retained  are,  in  fact,  PIC  codes,  discussed  in  a  subsequent  section.) 

The  WIRE  code  at  Lawrence  Livermore  National  Laboratory  is  a  good  example.®^  It  is  used 
routinely  to  compute  ETA  and  ATA  beam  transport  in  various  focusing  systems.  (The  ETA 
accelerator,  now  disassembled,  was  a  8  kA,  5  MeV  prototype  for  ATA.)  Among  the  elements 
considered  are  solenoidal  magnetic  lenses,  ion  channels,  grounded  wires  centered  on  the  axis,  and 
metal  foils.  For  the  centered  conducting  wire  the  radial  electric  field  is  fiixed  by  demanding  that 
the  potential  between  the  wire  and  the  drift  tube  be  zero,  which  requires  a  simple  integration  over 
the  disk  density  profile.  Fields  associated  with  focusing  foils  were  derived  by  Adler. A  typical 
1280  particle  run  with  numerous  focusing  elements  over  600  cm  takes  about  30  sec  on  a  CRAY-1 
computer. 

DISC,  a  similar  program  at  Sandia  National  Laboratories,  has  been  used  to  optimize  matching 
of  the  RADLAC  beam  onto  an  IFT  channel. The  experimental  configuration  modeled  is  shown 
in  Fig.  6.  A  20  kA,  16  MeV,  annular  electron  beam  enters  from  the  left  in  vacuum  along  a  17  kG 
axial  magnetic  field.  The  beam  was  born  at  an  immersed  cathode  and  so  has  significant  canonical 
angular  momentum;  it  is  in  the  slow  mode  of  rotation.  The  beam  then  enters  a  region  of  rapidly 
decreasing  magnetic  field,  where  it  spins  up  and  begins  to  expand.  The  expansion  is  arrested  by 
a  pair  of  focusing  foils,  and  the  beam  finally  is  captured  by  the  IFT  channel  at  the  right. 

Figure  7  is  a  sequence  of  three  particle  plots  illustrating  the  evolution  of  an  initially  off- 
center  beam.  Little  of  the  beam  current  is  lost  at  the  transition  despite  the  large  initial  offset. 
-Moreover,  the  beam  is  centered  by  the  IFT  channel  within  100  cm  of  transport,  although  its 
emittance  increases  substantially.  DISC  computations  such  as  this  consume  only  a  few  minutes 
of  minicomputer  time. 
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t  FIG.  6.  Experimental  geometry  for  DISC  code  simulation  of  RADLAC  beam  extraction  from 

magnetic  guide  field  and  capture  by  IFT  channel. 
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’.  Cross  sections  of  an  initially  offset  RADLAC  beam  at  (a)  0  cm,  (b)  40  cm,  and  (c)  100  cm 
after  extraction  from  a  17  kG  guide  field  into  an  IFT  channel,  computed  by  the  DISC 
code  for  the  geometry  depicted  in  Fig.  5.  The  outer  dashed  circle  is  the  drift  tube,  and 
the  inner  is  the  beam  RMS  radius  relative  to  the  drift  tube  axis. 
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Similar  single  disk  codes  have  been  employed  to  study  electron  beam  expansion  in  air  due 
to  scattering^’  The  Moliere’®  or  Keil-Zeitier-Zinn’'  formalism  should  be  utilized,  depenaing  on 
whether  many  or  few  scattering  events  occur  in  a  betatron  wavelength. 

III.  LINEAR  STABILITY  COMPUTATIONS 


A.  ^  Dispersion  Relation  Codes 

If  the  beam  equilibrium  is  not  too  complicated,  it  usually  is  possible  to  describe  the  growth 
of  linearly  unstable  modes  by  a  dispersion  relation  Fi^'.k,a].  where  is  the  (complex)  mode 
t'requency,  k  is  its  wavenumber,  and  a  the  set  of  parameters  describing  the  beam  and  accelerator. 
Peak  growth  rates  of  absolute  Instabilities  are  obtained  as  /m(_  )  for  the  zeroes  of  F.  mtiximized 
over  Peak  spatial  growth  rates  of  convective  instabilities  sometimes  are  given  by  the  temporal 
growth  rates  divided  by  the  corresponding  group  velocities,  d'^'idk.  again  maximized  over  k-  More 
often,  a  saddle  point  analysis  is  required.  The  convective  growth  of  the  resistive  wall  instability 
:r.  an  induction  accelerator  can  be  derived  by  this  method,  for  instance.’^ 

Complex  zeroes  of  locally  analytic  functions  are  easily  obtained  iteratively  by  .Muller  s 
method.’’^  In  essence,  the  algorithm  fits  a  quadratic  curve  to  three  values  of  the  function  and 
then  employs  the  curve  to  estimate  the  location  of  a  zero.  .\t  each  iteration  the  function  is  evalu¬ 
ated  at  the  estimated  zero  from  the  preceding  cycle,  and  a  new  estimate  made  based  on  this  and 
the  two  preceding  function  evaluations.  .A.  zero  typically  can  be  determined  to  a  relative  accuracy 
of  10'’  by  this  method  in  six  to  ten  cycles.  Subroutines  implementing  Muller’s  method  are  w'ideiy 
av-aiiable,  one  can  be  obtained  from  the  author.”^ 

.Although  Muller’s  method  can  be  applied  to  polynomial  functions,  the  Vieta  formula  solved  by 
•he  Newton-Raphson  method  is  more  efficient.®^  Moreover,  it  gives  all  the  zeroes  simultaneously. 
Many  implementations  are  available,  some  accommodating  repeated  zeroes.®^  There  are,  of  course, 
many  more  tecnniques  for  finding  zeroes  of  functions  than  the  two  mentioned  here. 

The  dispersion  function  need  not  be  known  analytically  to  apply  Muller's  method,  if  it  can 
oe  computed  numerically  with  sufficient  accuracy.  .An  example  is  the  GRADR  code,  which 
determines  the  normal  modes  of  cylindrically  symmetric,  radially  inhomogeneous,  particle  beams 
;r.  the  laminar  flow  approximation.®^  The  cylindrical  equilibrium  is  evaluated  from  six  nonlinear, 
coupled,  algebraic  and  first-order  differential  equations.  Four  follow  from  Majtwell's  equations 
oius  the  fluid  equations,  while  two  can  be  specified  arbitrarily  to  select  a  desired  equilibrium.'® 
The  two  constraint  equations  should,  of  course,  be  physically  realizable.  One  constraint  almost 
always  is  appropriate:  total  particle  energy,  kinetic  plus  potential,  must  be  constant  across  the 
oeam  and  equal  to  the  injection  energy.  .As  a  second  constraint,  typically  we  let  the  current  density 
profile  of  the  beam  be  specified  at  injection.  Conservation  of  canonical  angular  momentum  then 
determines  the  current  density  profile  within  the  accelerator  drift  tube.  The  resulting  system  of 
•'quations  IS  solved  iteratively.  Convergence  is  achieved  after  a  few  dozen  cycles,  unless  the  beam 

very  near  tne  space  charge  limit. 


Given  an  equilibrium  determined  in  this  or  any  other  way.  GR.\DR  then  solves  tne  ■ 
sponding  linearized  equations  to  obtain  eigenmodes  and  eigenvalues.  The  linear  equation  :'orm  .i 
fourth-order  Sturm-Llouvllle  system  In  radius.  The  equations  are  integrated  twice  from  t.ne  ax;.- 
to  the  outer  wall  with  distinct  arbitrary  choices  for  the  eigenfunction  values  at  the  axis,  .suocc: 
to  the  two  boundary  conditions  there.  .\  Fehlberg  fourth-fifth-order  Runge-Kutta  routine  " 
employed.  .Arbitrary  linear  combinations  of  the  two  solutions  then  are  used  to  evaluate  the  two 
boundary  conditions  at  the  wall,  and  the  determinant  for  the  unknown  coefficients  serves  as  t.ne 
dispersion  function  to  which  Afuller's  method  is  applied.  Solutions  can  be  obtained  at  the  rate 
one  per  second  on  a  CRAY-1  computer. 

GR.ADR  was  used  extensively  in  designing  the  R.ADLAC  accelerator. Typical  frequencies 
and  growth  rates  for  the  resistive  wall  instability  are  shown  in  Fig.  8.  Note  the  existence  r; 
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FIG.  8.  Predicted  (a)  frequencies  amd  (b)  growth  rates  of  the  m  =  1  transverse  resistive  wail 
instability  of  a  100  kA,  25  MeV  annular  electron  beam  propagating  along  a  20  kG  guide 
field  in  a  2.5  cm  radius  stainless  steel  drift  tube.  Straight  lines  in  (a)  are  the  edges  of 
cyclotron  and  space  charge  resonance  bands  due  to  the  variation  of  beam  parameters 
with  radius. 

branch  points  straight  lines  in  Fig.  8fa),  due  to  radial  inhomogeneity  in  the  beam  density.  High 
accuracy  must  be  demanded  of  the  differential  equation  solver  when  zeroes  lie  near  branch  points, 
as  is  the  case  here. 

GR.ADR  originally  was  written  to  investigate  collective  ion  acceleration.  In  many  collective 
acceleration  concepts  the  wave  phase  velocity  is  controlled  by  adiabatically  varying  some  external 
parameter  a.  such  as  the  guide  field  strength®*  or  the  drift  tube  radius.®®  The  corresponding 


rr.ange  in  wave  amplitude  can  oe  calculated  !rom  conser\'ation  ot'  wave  energy,  wn.cn  is  an  ootion 
n  the  code. 

D.  Linearized  PIC  Codes 

The  stabiiitv  ot  more  comipiicated  equilibria,  particuiariy  those  in  which  kinetic  effects  are 
important,  can  oe  established  with  linearized  PIC  codes.  In  such  models  equilibrium  fields  and 
particle  orbits  are  obtained  from  nonlinear  PIC  codes  or  by  other  means  le.g.,  the  ORBIT  code, 
.iiscussed  earner:.  Then,  perturbations  to  the  equilibrium  orbits  due  to  wave  fields  are  accumu¬ 
lated  aiong  the  unperturbed  orbits  and  used  as  source  terms  for  the  wave  fields  themselves.  This 
leap-frog  procedure  in  time  is  similar  to  that  of  standard  PIC  codes. 

The  linearized  particle  equations  of  motion  are  much  more  complicated  than  their  noniinear 
counterparts,  however,''  Formally, 

P,  -  E,  -  V-  ^  B,  -  r,  .  B.  -  A',  ■  c  [E.  -  Vh  ^  Bo)  <5,V  .  1C 


I'l  =  B,  ~  Py  P,  P-  (151 

-V,  (16) 

.\  savings  is  achieved,  therefore,  only  if  significantly  fewer  particles  are  needed,  and  the  effective 
dimensionality  of  the  problem  is  reduced. 

The  linearized  PIC  code  KMRAD,  for  example,  realizes  both  of  these  savings.’*  Cylindrically 
symmetric,  possibly  slowly  evolving,  beam  equilibria  are  determined  by  a  one-dimensional,  radi¬ 
ally  resolved,  noniinear  PIC  code  embedded  in  KMR.\D.  Three-dimensional  linearized  quantities 
aiso  are  resolved  radially  on  a  spatial  mesh  but  are  Fourier  decomposed  in  z  and  d\  one  (it,,  mi 
mode  IS  treated  at  a  time.  The  instability  growth  rate  of  the  fastest  unstable  wave  with  the  se- 
.ected  mode  numbers  is  determined  from  the  exponential  growth  of  field  energy  in  the  simulation. 
Field  and  current  radial  profiles  also  are  available.  Thus.  K.MR.\D  simulates  three-dimensional 
..nearized  dynamics  but  requires  computer  resources  comparable  to  those  of  a  one-dimensional 
rode.  Of  order  1000  particles  are  followed,  .\nswers  are  obtained  with  one  to  two  minutes  of 
CR.\Y-1  com,puter  time. 

KMR.\D  W2is  employed  recently  to  predict  two-stream  instability  growth  rates  in  IFT  channels 
for  recirculating  accelerators.’’  comparison  between  KMR.\D  kinetic  and  cold-beam  analytical 
results  for  a  cnannel  over-dense  by  a  factor  of  two  is  in  Fig.  9.  .\greement  with  growth  rates  from 
t.ne  l\’ORY  three-dimensional  PIC  code,  described  later,  is  good.  In  under-dense  channels  the 
eiectron-eiectron  instability  is  replaced  by  a  much  slower  electron-ion  instability.  The  code  aiso 
has  Deen  exercised  on  resistive  instabilities  in  much  higher  density  charmels. 
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FIG.  9.  Growth  rates  of  the  m  =  1  two-stream  instability  between  a  450  A,  35  MeV,  0.25  cm 
radius  electron  beam  and  electrons  of  an  1  cm  radius  IFT  channel  with  line  density 
initially  twice  that  of  the  beam. 


C.  Single-Component  Beam  Centroid  Codes 


The  previous  two  subsections  described  methods  of  obtaining  instability  growth  rates  with 
good  accuracy  at  fixed  points  in  an  accelerator.  Often,  it  is  more  desirable  to  obtain  the  total 
instability'  growth  with  less  accuracy  as  the  beam  propagates  the  length  of  the  accelerator.  This 
necessitates  a  simplified  description  of  the  beam  dynamics.  The  beam  breakup  instability  and 
most  other  m  =  1  modes  take  the  form  at  small  amplitudes  of  rigid  transverse  displacements  of 
the  beam.  Consequently,  the  beam  can  be  represented  as  a  string  of  rigid  disks  moving  forward  at 
a  specified  velocity  but  free  to  move  sideways.  Interaction  with  cavity  modes  and  other  sources  of 
instability  likewise  can  be  represented  simply  by  lumped  parameter  models.  Self-fields  are  treated 
in  the  long  w-avelength  approximation. 

The  beam  centroid  equation  of  motion  is 


d  d  d  2u 

_ .s _ '^  =  1 _  •  c _ c  ^  F 

dt  '  dt'‘  dt  '''  ' 


(1") 


wnere  f  =  .Vi  -  il'i  measures  small  transverse  displacements.  The  first  term  on  the  right  side 
)f  Eq.  1  I7'i  represents  the  magnetic  guide  field.  The  second  term  represents  the  net  image  force, 
eiectnc  less  magnetic,  for  a  cylindrical  drift  tube  of  radius  R.  The  final  term  includes  forces  from 
induction  modules,  etc.  Note  that  *?  increases  as  the  beam  accelerates. 

The  oscillatory  force  on  the  beam  due  to  m  =  1  cavity  modes  giving  rise  to  the  beam  breakup 
instability  is  given  by'^ 


T'h  -  —  —F  —  u/XF  =  I'ujr.  — 
dt  Q  dt  °  °  Q 


(18) 


A  separate  equation  is  required  for  each  induction  gap,  and  the  force  is  applied  impulsively 
there.  Z _  is  the  coupling  impedance,  Q  the  mode  quality  feictor,  and  cuq  the  mode  frequency.  The 
mage  displacement  instability  arises  from  interruption  of  beam  image  currents  at  the  acceleration 


gaps.  ' 


f  =  (19) 

For  narrow  gaps.  I  is  the  gap  width.  Otherwise,  /  must  be  computed  numerically  and  is  less  than 
the  gap  width. Finally,  finite  wail  conductivity  introduces  a  phase  lag  between  image  charges 
and  currents,  resulting  in  a  distributed  force, 


(20) 


derived  in  the  paraxial  approximation,  o  is  the  wall  conductivity,  and  g  is  a  geometrical  factor 
of  order  unity.  Evaluating  the  integral  is  relatively  time-consuming. 

Codes  employing  these  models  were  used  to  design  the  ATA^^  and  RADLAC“  accelerators. 
Figure  10  is  the  predicted  amplification  of  the  beam  breakup  instability  for  a  5  kA,  20  MeV',  19 
module  accelerator  design,  obtained  using  the  BALTIC  code.’^®  .Note  that  the  instability  wave 
packet  starts  at  the  beam  head  and  slowly  convects  back.  Growth  becomes  exponential  only  late 


;n  t;me.  This  calculation  required  nearly  a  minute  on  a  CRAV-1  computer.  The  code  wou:U  oe 
much  faster  with  the  resistive  wall  force.  Eq.  (20),  omitted. 

D.  Multi-Component  Beam  Centroid  Codes 

To  the  extent  that  focusing  forces  on  the  beam  particles  (including  self-forces)  are  not  linear 
with  radius,  particle  oscillation  frequencies  are  a  function  of  amplitude.  Phase-nux  damping 
of  transverse  instabilities  due  to  the  spread  in  particle  oscillation  frequencies  can  be  modeled 
within  the  context  of  beam  centroid  codes  by  treating  each  beam  slice  as  a  modest  number  of 
interpenetrating  rigid  disks.  Disks  are  assigned  effective  masses  and  charges  corresponding  to  the 
frequency  at  which  each  oscillates  and  to  the  portion  of  the  beam  each  represents,  respectively. 

Multi-component  beam  centroid  codes  are  used  to  simulate  electron  beam  hose  instabilities 
in  IFT  channels*’  and  in  air.  In  the  latter  context  they  are  described  at  length  in  the  chapter  by 
\V.  Fawley. 


IV.  NONLINEAR  PIC  COMPUTATIONS 


A.  PIC  Code  Fundamentals 

PIC  codes  determine  the  time  evolution  of  complex  beams  and  plasmas  by  computing  the 
dynamics  of  many  thousands  of  representative  particles  (electrons  and/or  ions)  moving  in  elec¬ 
tromagnetic  fields  externally  applied  or  produced  by  the  beams  and  plasmas  themselves.  Thus. 
PIC  codes  provide  the  most  fundamental  and  detailed  representation  possible  of  plasma  prob¬ 
lems.  In  effect,  they  solve  the  VHasov  equation.  Of  course,  this  precision  comes  at  the  cost  of 
substantial  computer  requirements,  and  for  this  reason  PIC  codes  should  be  employed  only  when 
simpler  numerical  or  analytical  techniques  described  previously  are  inadequate. 

The  electromagnetic  fields  are  defined  on  a  regular  mesh  in  one,  two,  or  three  dimensions, 
depending  on  the  symmetry  of  the  problem  to  be  solved.  The  mesh  can  be  in  rectangular, 
cylindrical,  or  other  desired  geometry.  .\t  each  time  step,  new  electric  and  magnetic  fields  are 
computed  by  advancing  the  finite  difference  approximations  to  Maxwell’s  equations, 
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using  currents  [J)  determined  from  the  plasma  particle  motion  on  the  previous  time  step.  \\- 
ternatively,  equations  for  the  scalar  and  vector  potentials  can  be  solved.  Boundary  conditions 
are  required  to  define  spatial  derivatives  at  the  mesh  edges.  Wave  reflecting  (i.e.,  metallic)  or 
periodic  boundary  conditions  are  common  choices.  More  complicated  boundaries  allow  electro¬ 
magnetic  waves  to  be  launched  into  the  computational  region  or  to  leave  it.  as  described  in  the 
next  subsection. 


Particle  momenta  [P  =  ^V),  positions  (A’),  and  energy  ("rj  are  then  advanced  using  the 
relativistic  equations  of  motion  with  the  newly  computed  fields. 

P  =  E-VxB  i23) 


A'  =  V' 


\  t  ** 

=  (p=-‘)  • 


(Replace  by  1  for  nonrelativistic  problems.)  The  fields  appearing  in  Eq.  (23)  are  those  at  the 
particle  location,  obtained  from  the  fields  at  nearby  mesh  points  by  (typically,  linear)  interpo¬ 
lation.  When  a  particle  leaves  the  computational  region,  it  is  destroyed  or  it  is  returned  to  the 
mesh  by  some  prescribed  procedure  (e.g.  reflected).  By  the  same  token,  particles  can  be  injected 
from  boundaries,  a  feature  particularly  useful  in  particle  beam  simulations.  After  a  particle's  new 
position  and  momenta  have  been  determined,  its  contribution  to  the  plasma  currents  is  obtained 
by  interpolating  V’  =  P  to  nearby  mesh  points. 

This  cycle  of  advancing  fields  based  on  particle  currents  and  then  advancing  particles  based 
on  the  new  fields  is  repeated  hundreds  of  times  in  a  typical  simulation.  The  time  step  is  set  by  the 
smallest  time  scale  in  the  calculation,  which  may  be  the  plasma  oscillation  period,  the  electron 
cyclotron  period,  or  the  Courant  time  (of  order  the  time  for  a  light  wave  to  cross  a  cell  in  the 
mesh).  Progress  has  been  made  in  the  last  few  years  at  surmounting  the  Courant  limit,  which  is 
numerical  rather  than  physical  in  character. Cell  dimensions  must  be  small  compared  to  spatial 
scales  of  interest. 

PIC  code  running  times  and  memory  requirements  are  highly  problem  dependent.  CPU  times 
on  a  CRAY-1  computer  typically  range  between  15  minutes  and  4  hours,  although  20  hour  runs 
are  not  unheard  of.  Corresponding  central  memory  needs  vary  between  2- 10^  and  4- 10*’  words. 
.\t  least  two  fast,  large  capacity  disks  or  their  equivalent  also  are  needed.  Historically,  the  physics 
problems  attempted  with  PIC  codes  have  expanded  to  consume  the  maximum  resources  available 
in  each  generation  of  computers. 

PIC  codes  usually  have  extensive  graphics  output  capabilities  and  operating-system  interfaces. 
.\dvances  in  PIC  technique  are  described  in  several  books**  **  as  well  as  in  the  chapters  here  by 
.\.  Langdon  and  A.  Mankowsky. 


B.  Standard  PIC  Codes 

Most  any  multidimensional  PIC  code  treating  both  electric  and  magnetic  fields  can,  with 
minor  modifications,  be  applied  to  beam  transport.  Beam  simulations  typically  are  performed 
m  cylindrical  geometry.  Relativistic  particle  dynamics,  if  required  and  not  already  available,  are 
easy  to  implement.  Efficiency  is  an  important  consideration:  beam  simulations  tend  to  be  lengthy. 


Boundary  conditions  properly  modeling  the  often  complex  accelerator  structures  are  partic¬ 
ularly  important.  Figure  11.  for  Instance,  is  taken  from  an  IVORY  code  simulation  of  an  image 
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FIG.  11.  .Axisymmetric  drift  tube  geometry  from  FV'ORY  simulation  of  image  displacement  insta¬ 
bility  with  23  kA,  2  MeV',  annular  electron  beam  in  a  10  kG  guide  field.  Beam  dynamics 
are  nonaxisymmetric. 

displacement  instability  experiment. Shaded  regions  are  metal  disks:  thick  horizontal  lines  pro- 
.ecting  from  the  disks  are  metal  foils.  Conductors  are  introduced  on  the  computational  mesh 
simply  by  setting  tangential  electric  fields  to  zero  on  their  boundaries.  (The  other  field  com¬ 
ponents  need  not  be  adjusted  on  a  properly  staggered  mesh.)  A  tightly  spaced  helical  coil  can 
be  specified  by  setting  the  field  component  along  the  winding  to  zero  while  leaving  the  other 
components  alone. Specifically,  for  a  helix  of  pitch  angle  v,  set  cos  v  i?,  sin  v  =  0  without 
modifying  £■,  sinv  -  £,  cosv.  Dielectrics  are  treated  by  multiplying  E  in  Eq.  (21)  by  the  dielec¬ 
tric  constant  «.  Similarly,  bulk  resistivity  is  included  by  adding  the  term  oE  to  the  left  side  of 
Eq.  (21).  In  so  doing,  one  must  be  careful  not  to  exceed  the  Courant  limit  by  allowing  cAt  to 
become  large,  unless  integrating  factors  are  used. 

If  the  simulation  mesh  is  not  large  enough  to  accommodate  the  entire  accelerator,  then  com¬ 
putational  boundaries  through  which  electromagnetic  waves  are  able  to  pass  without  nonphysical 
reflections  are  essential.  The  standard  wave-transmitting  boundary  condition  is 

dE  dE  dEo  ,  , 

=  '“I 

where  E  is  the  tangential  electric  field  at  the  boundary  and  x  is  the  normal  coordinate,  pointing 
inward.  Any  electromagnetic  wave  approaching  the  boundary  with  normal  ph2ise  velocity  v  leaves 
the  mesh  without  reflection,  as  can  be  verified  by  substituting  arbitrary  E{t  x/v]  into  Eq.  (26). 
Equation  (26)  also  launches  waves  Eoit  —  x!v)  into  the  computational  mesh,  if  desired.  When 
the  primary  source  of  electromagnetic  waves  is  short  wavelength  beam  fluctuations,  setting  v  to 
the  beam  velocity  minimizes  spurious  noise. 

If  waves  with  a  spread  in  phase  velocities  strike  the  boundary,  at  least  some  of  them  are 
reflected.  The  reflected  wave  relative  amplitude  is  (u*  —  u)/(u'  -t-  v),  where  v'  is  the  wave  phase 
velocity  normal  to  the  boundary.  A  more  general  boundary  condition  able  to  transmit  waves  at 
essentially  all  phase  velocities  was  developed  for  laser-plasma  intereiction  studies.**  It  is  relatively 


complex  and  will  not  be  discussed  here.  A  second  alternative  for  wave  absorption  is  a  tinicK, 
graded  resistance  placed  at  the  boundary,  although  a  large  portion  of  the  computational  mesn 
can  be  lost  in  this  way. 

-\n  example  of  multidimensional  simulations  of  accelerators  is  a  recent  instability  study  for  high 
current  modified  betatrons,  performed  with  the  IVORY  code.*®  fV'ORY  is  a  three-dimensional, 
relativistic,  electromagnetic  PIC  code  in  which  fields  are  Fourier  decomposed  in  the  azimuthal  di¬ 
rection  rather  than  finite-differenced.*^  This  approach  is  particularly  efficient  for  problems  which 
are  approximately  but  not  exactly  axisymmetric,  such  as  transverse  oscillations  on  an  axisymmet- 
ric  beam  or  toroidal  oscillations  on  an  electron  ring.  Applied  to  the  modified  betatron,  IPROP 
can  follow  instabilities  in  the  recirculating  beam  for  100  m  in  a  few  hours  of  CRAY-1  computer 
time.  \  typical  negative  mass  instability  result  is  shown  in  Fig.  3. 

The  application  of  PIC  codes  to  ion  beam  equilibrium  and  stability  is  discussed  in  the  chap¬ 
ter  by  W.  Herrmannsfeldt.  Typically,  a  two-dimensional  code  is  used  to  follow  the  transverse 
dynamics  of  a  single  beam  slice  as  it  travels  the  length  of  the  ciccelerator. 

C.  Darwin  Model  PIC  Codes 

.\s  noted  in  Sec.  I\’.  .  the  time  step  in  electromagnetic  PIC  codes  is  limited  by  the  Courant 

condition.  This  limitation  can  be  circumvented  by  implicit  methods,  some  of  which  are  described 
in  the  chapter  by  A.  Langdon.  Alternatively,  electromagnetic  waves  can  be  eliminated  entirely,  if 
they  are  not  necessary  to  the  problem  to  be  solved.  Electrostatic  calculations  are  entirely  adequate 
for  many  ion  beam  calculations,  for  instance.  The  Darwin  model  can  be  applied  efficiently,  when 
inductive  effects  are  important  but  wave  effects  are  not. 

Darwin's  approximation  to  .Maxwell’s  equations  is  derived  formally  by  expanding  the  particle- 
field  interaction  Lagrangian  in  powers  of  V/c,  where  V  is  the  particle  velocity,  and  retaining 
terms  through  second  order.  The  field  equations  following  from  the  approximate  Lagrangian  can 
be  written  in  several  forms,  including** 

El  =  -Vd)  V*d)  =  -47rp  (27) 

(28) 

ot 

d 

^  <  B  =  47rJ  ^  —  El  (291 

dt  ^ 

The  subscripts  and  “f"  designate  longitudinal  and  transverse  components  of  the  fields  and 
currents;  V  ■  E,  =0.  Poisson’s  equation  must  be  solved  to  obtain  J,  and,  of  course,  <?>.**  Nonethe¬ 
less.  the  added  cost  of  solving  Eqs.  (27)-(29)  instead  of  Eqs.  (21)-(22)  is  fully  justified  by  the 
larger  time  step  permitted.  .Numerical  noise  is  less  also. 

The  Darwin  model  is  implemented  in  the  two-dimensional  DPC  code,  used  to  refine  beam 
injection  and  transport  systems  in  the  ,AT.\.^  .\n  interesting  issue  explored  with  DPC  is  acceler- 


ation  of  electrons  from  the  IFT  channel  used  to  guide  and  stabilize  the  primary  beam.  Figure  1 
from  this  study,  shows  the  accumulated  channel  electron  current  and  energy  at  induction  mod'. 
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FIG.  12.  Channel  electron  current  (a)  and  energy  (b)  at  AT.\  gap  50.  predicted  by  the  DPC  code 
under  certain  circumstances. 


D.  PIC  Codes  with  CoMoving  Meshes 

When  performing  simulations  of  beams  propagating  in  plasma  channels,  it  often  is  convenient 
to  use  a  computational  mesh  moving  with  the  beam.  Conceptually,  this  is  accomplished  by  a 
Lorentz  or  a  Galileam  transformation  to  the  beam  frame.  After  a  Lorentz  transformation,  PIC 
codes  of  the  sort  discussed  above  can  be  employed.  A  Galilean  transformation  is  preferable, 
however,  when  a  resistive  medium  is  present  or  when  the  range  of  time  scales  can  be  comp'^essed. 

The  electromagnetic  field  equations  in  a  coordinate  frame  moving  axially  at  velocity  v  are 
given  by®^ 


n  -  -  ''si  ^  [I  -  -  '''si  '=■  - 


5  [s  - 1'  -  ''ll  -  "'ll  =  1^- 
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L.gnt  wave  axial  characteristics  in  this  frame  are  represented  schematically  in  Fig.  13.  Note  that 
the  forward-going  characteristic  crosses  a  constant-Z  grid  line  before  it  crosses  a  constant-T  grid 


line,  if  Az  <  !c  —  i')At.  Information  at  (i  —  l)A2  influences  the  flelds  at  lA^  at  the  same  time 
.evei  in  this  case,  and  the  Courant  limit  is  relaxed  somewhat.  If,  in  addition,  the  electromagnetic 


coupling  among  radial  positions  at  fixed  Z  and  T  is  treated  implicitly,  the  Courant  limit  becomes 
Af  <  A;  ic  -  v).  The  simulation  time  step  typically  is  constrained  by  particle  dynamics  under 


tnese  circumstances. 


This  algorithm  is  implemented  in  a  version  of  the  IVORY  code  known  as  IPROP.''‘  IPROP 
IS  applied  routinely  to  studies  of  relativistic  electron  beam  propagation  in  low  and  high  den¬ 
sity  plasma  channels.*®  .A.  lumped-parameter  plasma  conductivity  package  is  available  to  treat 
weakly  ionized  gases  in  an  Ohm’s  law  approximation.  Particle  scattering  and  bremsstrahlung 
routines  also  are  provided. 

significant  savings  in  computer  memory  and  a  moderate  savings  in  computer  time  can 
be  achieved  for  simulations  of  ultrarelativistic  electron  beams  by  employing  the  frozen  field 
approximation.^’  The  frozen  field  equations  are  obtained  from  Eqs.  (30)-(35)  by  setting  v  =  c 
and  dropping  time  derivatives.  Fields  are  solved  at  each  particle  time  step  by  integrating  the 
equations  in  Z  from  the  head  of  the  beam  to  the  tail.  Application  of  the  frozen  field  equations  to 
simulations  of  beam  propagation  in  air  is  reviewed  in  the  chapter  by  W.  Fawley.  The  frozen  field 
approximation  is  useful  in  analytical  work  as  well,  e.g.,  for  the  plasma  wakefield  accelerator.®  Be 
aware,  however,  that  the  approximation  fails  for  beams  passing  through  metal  foils  or  apertures. 
Like  the  Dairwin  approximation,  the  frozen  field  approximation  does  not  treat  electromagnetic 
radiation.^ 
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